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Abstract The hard sphere repulsion among ions can be considered in the Poisson-Nernst-Planck 
(PNP) equations by combining the fundamental measure theory (FMT). To reduce the nonlocal com¬ 
putational complexity in 3D simulation of biological systems, a local approximation of FMT is derived, 
which forms a local hard sphere PNP (LHSPNP) model. In the derivation, the excess chemical poten¬ 
tial from hard sphere repulsion is obtained with the FMT and has six integration components. For the 
integrands and weighted densities in each component, Taylor expansions are performed and the lowest 
order approximations are taken, which result in the hnal local hard sphere (LHS) excess chemical 
potential with four components. By plugging the LHS excess chemical potential into the ionic flux ex¬ 
pression in the Nernst-Planck equation, the three dimensional LHSPNP is obtained. It is interestingly 
found that the essential part of fre e energy te rm of the previous size modified model(|B orukho v et all 
Phys. Rev. Lett. 1997, 79, 435-438; Kilic_eT_d Phys. Rev. E 20 07, 75, 021502: ILu and ZhorJ Biovhvs. 
J. 2011, 100, 2475-2485: iLiu and Eisenberd J. Chem. Phys. 2014, I 4 I, 22D532) has a very similar 
form to one term of the LHS model, but LHSPNP has more additional terms accounting for size effects. 
Equation of state for one component homogeneous fluid is studied for the local hard sphere approxi¬ 
mation of FMT and is proved to be exact for the first two virial coefficients, while the previous size 
modified model only presents the first virial coefficient accurately. To investigate the effects of LHS 
model and the competitions among different counterion species, numerical experiments are performed 
for the traditional PNP model, the LHSPNP model, the previous size modified PNP (SMPNP) model 
and the Monte Carlo simulation. It’s observed that in steady state the LHSPNP results are quite 
different from the PNP results, but are close to the SMPNP results under a wide range of boundary 
conditions. Besides, in both LHSPNP and SMPNP models the stratification of one counterion species 
can be observed under certain bulk concentrations. 


Keywords Hard sphere repulsion • Three dimensional fundamental measure theory • Poisson-Nernst- 
Planck equations ■ LHSPNP model • Equation of state 


1 Introduction 

In systems of charged particles, the traditional continuum Poisson-Nernst-Planck (PNP) model has 
been widely used for modeling and simulating a variety of interesting phenomena, such as ion distribu- 
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tions around charged particles [1, [ 13 , , and flow of ions across channels d H, [llllO, HHli, IH and 
axons [i^. A defect in the PNP model is the point charge treatment of ions which in fact have finite vol¬ 
umes. This may lead to unphysical concentration values in the neighborhood of biomolecules [H ll^. l4l| 
and hence can not properly capture the size effects which may play important role in ion permeation 
through ion channels. 


A great deal of efforts have been devoted to remedy this disadvantage and improve the PNP 
model to provide reliable predictions through including ionic size effects and ion-ion correlations. In 
equilibrium and uniform ionic size, Andelman et al. have incorporated the ionic size effects into the 
traditional continuum model, and provided the explicit expressions of ion concentrations and the 
uniform size modified Poisson-Boltzmann (SMPB) equationjl, Q. This is derived by introducing an 
additional solvent entropy term, representing the unfavorable energy modeling the overpacking or 
crowding of ions and solvent molecules, into the free energy form. Fenley’s group has studied this 
SMPB model explicitly through comparisons with the original Poisson-Boltzmann (PB) model, and 
investigations of the sensitivities to parameterization |lM Il3 .1^ . By using a rigorous lattice gas method, 
Chu et al. extended this uniform SMPB model to make it work for two different ion sizes and studied 
ion binding to DNA duplexes which showed improved agreements with experimental data@. In the 
case of three or more different sizes, the rigorous lattice gas method is hard to apply to give an explicit 
SMPB equation, but an alternative size modified Poisson-Nernst-Planck (SMPNP) model is derived 
in our previous work which can be considered as an implicit SMPB model in equilibrium ^ . Li 
et al. have provided some mathematical analysis on the energy functional for the SMPB model 3l|. In 
recent work by Liu and Eisenberg [33 - I^ . they considered the size effects from all spherical ions, water 
molecules (treated as polarizable spheres) and interstitial voids, and derived an entropy form of those 
spheres of arbitrary diameter and voids. To further account for the correlation effect of ions, Liu and 
Eisenberg’s work also included a high-order derivative term of potential in free energy functional derived 
by SantangelojUj which has similar form to Cahn-Hilliard concentration-gradient expansions [M 1^. 
The resulted fourth order Poisson-Fermi equation [H, |35| (for equilibrium condition) and Poisson- 
Nernst-Planck-Fermi equations (for nonequilibrium condition) can account for the steric effect of 
ions and water molecules, the correlation effect of crowded ions, the screening effect of polar water, 
as well as the charge/space competition effect of ions and water molecules. These models were shown 
to be able to capture critical phenomena such as ion binding, blocking, and selective permeation of 
L-type calcium channels 


Apart from the above mean field modified models, coupling of density functional theory (DFT) and 
PNP model has also been used to account for ion-ion interactions in solutions. The major description of 
DFT lies in the construction of the excess Helmholtz free energy [l^. Rosenfeld derived the fundamental 
measure theory (FMT) to express the excess Helmholtz free energy for hard sphere mixtures [45^. which 
contains four scalar and two vector weight functions, and has been widely used in many literatures [TTl 
[T3 . [T3 . [TtI [s^, HI, m, HO, [^ l46i - l^ . m, IsqI - I^ 1131 • Not long after that, he extended the 

original FMT, that works for three dimensional inhomogeneous fluids, to predict a freezing transition 
of the hard-sphere fluid into a solid [Slj, [53|. At the same time, lots of contributions have been made to 
improve the FMT. Kierlik and Rosinberg have proposed a simplified version of FMT, which requires 
only four scalar weight functions [1^. Roth et al. and Wu et al. have independently, published the white- 
bear versionor the modified FMT[63] using the Mansoori-Carnahan-Starling-Leland bulk equation 
of state to make simulation results more accurate. For inhomogeneous fluids of nonspherical hard 
particles, Mecke et al. have derived a fundamental measure theory using the Gauss-Bonnet theorem[l3|. 
In numerical calculation, due to the complexity of FMT, calculations are mostly reduced to ID cases 
using geometric symmetry of the simulation systems[l7Ll25ll53l.l54ll64l|. Only a few literatures provided 
3D results [T3 ,[i3)[30j [53 ■ When FMT is combined with the PNP model to include hard sphere repulsion 
for non-equilibrium transport studies, the final integro-differential equations are more complex and 
computationally expensive, especially in the situation with biomolecular systems. Liu et al. have made 
some mathematical analysis on the PNP system for ion flow with density functional theory for hard 

S here interaction in ID situatio n [s^. So far, most studies are restrained to ID caseflj, [13 HU 
, [ 3 ^. Recently, Gillespie et al.[28l| and Meng et al.ld^ have reported simulations of 3D systems. 
The main difficulty is the calculation of the nonlocal components originated from the Helmholtz free 
energy. 


In this work, we propose a 3D local hard sphere PNP (LHSPNP) model for real biomolecular systems 
in ionic solutions. On one hand, 3D simulation can count for the irregular boundary of biomolecule 
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which plays an essential role during biological processes, such as ionic flow across channel, protein 
modification or interaction with a protein or substrate molecule and cell signaling . These shape 

information is hard to be captured in ID case resulted from the symmetric boundary simplifications. 
On the other hand, the local hard sphere (LHS) model largely simplifies the numerical calculations 
while it still partially maintains the effects of hard sphere repulsion. 

We improve the PNP model by considering the hard sphere repulsion from FMT. The excess 
Helmholtz free energy of FMT is employed to generate the excess hard sphere chemical potential 
which is completely ignored in the PNP model. Different from the ideal chemical potential, this excess 
component of a certain ion at a given point is determined by an integration about all ion concentrations 
in a region around the given point, rather than the certain ion concentration at the given point. This 
integration is in the form of convolution. Gillespie et al. have used fast Fourier transform to deal with 
the convolutionf^ . while Meng et al. have used the definition of Dirac delta function and change of 
variables to transform these 3D integrals into 2D integrals on spheres and remove the singularity in the 
integrands [d^. Though both algorithms can solve the integro-differential equations, they cost lots of 
computer memory and time during calculation. We aim to construct an excess chemical potential in a 
point to point way like the ideal component for 3D simulations of ionic solutions. Liu et al. have derived 
a LHS excess chemical potential in ID case and made theoretical analysis on the model problem based 
on geometric singular perturbation theory [^. In our work, we concentrate on the more complex 3D 
case and simplify the integration using an expansion of the integrand under small ionic diameters to 
obtain the final LHSPNP model. It’s found that the 3D LHSPNP model is exactly the same as Liu’s 
when reduced to ID case. 

We examine the effects of the LHS model both from a theoretical and a numerical point of view. 
For the case of one component homogeneous fluid, the virial coefficients are investigated through 
the equations of state from the FMT, our local approximation and the size modified model [1, [l^, [13] 
studied in former work. For FMT, the equation of state is known to be the Percus-Yevick compressibility 
equation, equivalent to scaled particle theory, which provides the first three virial coefficients exactly. 
Frydel and Levin have shown that the size modified model only predicts the first virial coefficient 
exactly [Till. Based on these discussions and the relations between bulk grand canonical potentials from 
the micro and macro views, we derive the virial coefficients for our local approximation. We find that 
the LHS model provides the first two virial coefficients exactly, which performs better than the size 
modified model. Numerically, to investigate the effects of the hard sphere repulsion from LHS model, 
we make numerical comparisons for a spherical cavity case among the PNP, LHSPNP, SMPNP models 
and the Monte Carlo (MC) simulation. For the four parts of the local excess chemical potential in 
LHSPNP, we notice that the first term is, to some extent, similar to the local excess chemical potential 
of the SMPNP model, though LHSPNP and SMPNP are based on two different theoretical frameworks. 
We consider two counterion species in the solution, which helps us to understand competitions between 
them under the PNP, LHSPNP and SMPNP models. 

The rest of the paper is organized as follows. The method section offers a detailed description about 
the derivation of the LHSPNP model, its equation of state in one component uniform fluid, and the 
numerical method to solve the LHSPNP equations. The result and discussion section provides numerical 
results of a spherical cavity in various bulk concentrations with the PNP, LHSPNP, SMPNP models 
and MC simulation and also some discussions about the observed phenomena. Finally, conclusions are 
summarized in the conclusion section. 


2 Method 

2.1 Local Hard Sphere Poisson-Nernst-Planck Model 

In addition to the ideal chemical potential considered in the original Poisson-Nernst-Planck (PNP) 
equations, an excess chemical potential arising from hard sphere repulsion is incorporated to the 
equations to make the model more accurate. Using the constitutive relations about the flux and the 
electrochemical potential, we have the following expression of ion flux 

Ji = -midV^ii = - 


CiVni, 


( 1 ) 
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where Ji is the flux, rrii and Ci are the mobility and concentration for the zth ion species respectively, 
fj.i is the chemical potential, Di is the diffusion coefficient, ks and T are the Boltzmann constant and 
absolute temperature. The chemical potential /ii is composed of two parts, the ideal part fj,"^ and the 
excess part /i®^. The mass and current conservation law leads to the following equation 

^ = -V • J. = V • = V • , (2) 

where = qicj) + InA^a, (j) is the electrostatic potential, /3 = Qi and Ai are the charge 
amount and de Broglie wave length of the zth ion species respectively. Coupling Eq.[2]with the Poisson 
equation, we obtain the modihed PNP at steady state 

V-A(Vc, + ^ftQV(^ + /3c,VMr) = 0, ^ = (3) 

K N 

-V ■ e\7(j)-'^c^qi = '^QiSi, (4) 

2=1 2=1 


where K is the number of ion species, and Qi is the charge amount of the ith atom in the biomolecule 
that contains N fixed point charges. Specifically, it is worth noting that in equilibrium state of zero 
flux, the above modified PNP equations become a modified PB model as described by the following 
equation; 


— V ■ e{x)\/(j){x) 




(5) 


where is the bulk excess chemical potential defined under bulk concentrations. This generalized 
PB equation incorporating the hard sphere repulsion can be solved by combining with the following 
local expression of ^f^{x) from FMT. 

According to Rosenfeld’s fundamental measure theory (FMT) [4^, the excess Helmholtz free energy 
due to hard sphere repulsion can be expressed as 


^ ec 


{cz(a^)} 

{na(a;)} 


= /3 


-1 


dx^ 


{na(a;)} 


= -no ln(l - ns) + 


1 


1 - no 


(nin2 - nyi • nv2) 


n2 


- 3nv2 ■ nv2), 


,{x) = 


247r(l -n3)2 
J Ci{x')uj'f\x — x')dx', 


( 6 ) 

(7) 

( 8 ) 


where na{x) is the weighted density, uj^°‘\x) is the characteristic (weight) function for a = 0,1,2, 3, VI, V2. 
The indexes V 1 and V 2 are used to represent the vector terms while the other four represent the scalar 
terms. For a three-dimensional hard sphere particle of radius Ri, these functions are defined as: 


Ujf\x) = e[\x\- Ri), (9) 

ujf\x) = 5{\x\-Ri), (10) 

uj^^\x) = uj^'^\x)/AttR„ ( 11 ) 

u}f^\x) = Jf'{x)/^'KR\, (12) 

= |^^(|a;| - Ri). (13) 

‘^i^'^Hx) = u}’i^^\x)/4:TtR,. (14) 

9{x) is the unit step function defined by 


0 {x) 


0 if X > 0, 
1 if X < 0, 


(15) 
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and S{x) is the Dirac delta function. 

From above definitions, it’s not difficult to obtain the following expression of the excess chemical 
potential 


where 


Mr( 


Sci 


= P 


-1 


dx 


^ drir. 




d<P 

dno 

d<P 

dni 

dn2 

d<d> 

dns 

d<P 

dnvi 

d4> 

dnv2 


-ln(l - ns), 

n2 

1 - ns’ 

, 1/2 ^ 

i- a~T\ -^(^^2 - '^V2 ■ nv2), 

I-ns 87r(l 

”0,1/ ^ , n2 

+ 1241 -ns)^ 

x(n2 - 3nv2 • nv2), 
nv2 
1 - ns ’ 

nvi 1 n2nv2 
1 — 7/3 477 (1 — nss)^ 


(16) 

(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 

( 22 ) 


As shown in Eq. 1161 the excess chemical potential resulted from hard sphere repulsion is an integra¬ 
tion defined on a certain region, which requires plenty of calculation. By taking corresponding Taylor 
expansions of the excess chemical potential, we can deduce a local formula for the excess chemical po¬ 
tential as follows. From Eq. [1^1 the excess chemical potential of the ith ion species can be decomposed 
into six components. Denote them by {//“(a;) = j5~^ J [{n-y(£c)}] w-“^(£c — x')dx'} and we can get: 



= - In (1 - E ^^R^jix)) + 0 {R^), 

(23) 


RiY\-AT:RpCj{x) 

1-Ej s^RjCjix) 

(24) 

/3m?(®) 

47ri?f JO - i?,c, (a:) 

= V/ 

1 - Ej s^RjCjix) 

(25) 


4 RiYZi ^i{x) A 

= -7: +OiR% 

3 1-Ej|7ri7jcj(a;) 

(26) 


= 0{R^), 

(27) 


= 0{Ry). 

(28) 


Here, we give the explicit derivation of Eq. [23] as an example. Others can be obtained in the similar 
way. First, consider the expansion of ns{x) in the definition of P^^(x)\ 

ns{x) = ^^ J Ci{x')uj^^\x — x')dx'= J Ci{x — x')ujf‘\x')dx' 

i i 

= / Ci{x — x')dx' 


'\x'\<R- 


= E / 

^ J\x'\<Ri 


Ci{x) — Vci(x) ■ x' -I- 0{B?) 


dx' 


= E o^-^?c,(a;) + 0{R^). 


(29) 
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Then, substituting above equation for n^{x) in the expansion of /3/i°(a;), we have: 

= J ^-{x')ujf\x - x')dx'= J ^-{x - x')u}f\x')dx' 

= - f \n{l- mix - x'))^^^d{\x'\ - Ri)dx' 

•J 

1-217 p17 poo 2 

= — d4> j sinOdO / \tl [l — n-^ix — rv)^ - ^S{r — Ri)dr 

Jo Jo Jo 47ri?j 

o'Z'K oTZ 

= -/ d<f> \n (l — nsix — Rii')) sin 9d9 

47r Jo Jo 


'0 ^0 
•27r p-K 


In (l - 713 ( 0 ;)) + • i-RiJ^) + 0 (||V^n 3 (a;)||i?^) 


477 70 ^^Jo ' 1 - 773 ( 3 ^) 

= - In (1 - 773 ( 0 ;)) + 0 {R‘^) 

= - In (1 - ^ |7ri?|c,(o;) + 0 {R^)) + OiR^) 


sin 9d9 


- In (1 - ^ -TTR^jCjix)) + OiR^), 


(30) 


The other terms can be obtained similarly, as shown in Eqs. Take the lowest order approxi¬ 

mation of each term and the local excess chemical potential in three dimensions is given by: 


/3/if^^(o;)=/3^^“(o;) 


-ln(l-^-7ri?|cj(o;)) 


Ri J2j ^TTR'jCj^x) 
1 - Ej |7ri?|cj (o;) 


-f 


47ri?2 . Rjcj{x) 


1 - Ej hRJjCj{x) 


4 R^iY.jCj{x) 

3^1 - ^TTR^-Cjix)' 


(31) 


It’s notable that the excess chemical potential expressed by Eq. [21] at a given point is determined by 
the concentration values at this point. This local expression is much simpler than the nonlocal one in 
numerical calculation. For one dimensional situation, a local hard sphere potential is also proposed by 
Liu et al.[^ to investigate ion flow through channels and shows great improvements. 

The final modified Nernst-Planck (NP) equations of the local hard sphere PNP (LHSPNP) can be 
obtained by replacing the /r™ in Eq. |2|by in Eg. HTH 

V • Di{x) (Vci{x) + PqiCi{x)\7(j){x) + l3ci{x)V (x)'^ =0, i = ,K. (32) 

In the size modified PNP (SMPNP) equations proposed in our former work [13,1431, the excess 
chemical potential is expressed by iif^{x) = —j3~^ki In (l — Y) ■ aJcj{x )\, where ki = at and oq are 
the diameters of ion and water molecule, respectively. This is quite similar to the first term of the local 
hard sphere excess chemical potential ij-^^^{x) in Eq. |2I1 From this point of view, the LHSPNP can 
capture the same ionic size effects contained in the SMPNP model. Furthermore, it’s found that all the 
four terms of Eq. |3T] are positive and contribute unfavorable energies for all possible concentrations, 
which indicates the later three terms have the similar influence as the first term and strengthen the 
hard sphere repulsion effect. Compared with the excess chemical potential in SMPNP model, these 
extra three terms in fif^^{x) can be regarded as supplements to the size effects. 
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2.2 Bulk Fluid Equation of State 

In uniform fluid with constant density, the equation of state for FMT is known to be the Percus- 
Yevick compressibility equation, which is the same as the scaled particle theory [T^ li^, [ 53 . The four 
scalar weighted densities {na}{a = 0,1,2,3) can be reduced to where 

^( 2 ) = 47rJ))- = ^^CiRi and while the vector weighted densities nyi and nv 2 

vanish. The bulk grand canonical potential is given by 

12b = + n‘'-J2 (33) 

i 

where l3~^<Pb is the bulk excess energy density defined in Eq. 0 is the bulk ideal energy density 
defined by /3“^ Ci(lnciAf — 1) and V is the system volume [l(ll|. On the other hand, 12b can also be 
determined from the thermodynamic relation 


12b = -PV, 


(34) 


where P is the pressure [l^- From Eg. IHHland Eq. [MJ the pressure is expressed as: 


i 


(35) 


According to the equilibrium condition, the chemical potential at bulk state is where / is the 
energy density composed of the excess and ideal parts in the EMT. This leads to the general expression 
of pressure in the following equation: 


i 



(36) 


In studying the equation of state, we consider one component bulk fluid with constant density p. 
Under this assumption, the bulk excess energy density from Eq. 0 is 


^6 


—pln(l — ??) + 


3??P 

1 -?7 


2 (l-r,)2’ 


(37) 


where ry is the packing fraction defined by 77 = and a is the diameter of ions in solution. Thus we 
can get the energy density, and the equation of state 


P 


1 + 77 + 77^ 
(1-77)3 


i+E 


37^ F 37 T 2 
2 


1 + 477 + 1077^ + 1977^ +_ 


(38) 


Based on tabulated values of the first eight virial coefficients [T^, the equation of state for homogeneous 
hard sphere fluid is given by 


^ = 1 + 477 + lOry^ + 18.355773 + 28.22577^ + 39.7477® + 53.577® + 70 . 877 ^ + .... (39) 

P 

It’s apparent that the first three virial coefficients are exact from FMT in Eq. |3H1 

In our local hard sphere approximation, we replace the weighted densities with a local expression 
by ignoring the Taylor expansion terms of order 0{a"^){m ^ 4), and use it in the derivation of excess 
chemical potentials. With this consideration, the bulk excess energy density in local hard sphere model 
is given by 


^LHS 


-pln(I - 77 ) + 


^VP 

1 - 77 ' 


(40) 








Compared with Eq. [321 this equation does not take into account the last term of order 0[rf) = 0[a^). 
Following this definition, we can get the bulk excess chemical potential after ignoring the terms of 
order 0(cr™)(m ^ 4) 




d^LHS 

dp 


= - ln(l -T]) + 

1—7] 


(41) 


This is in accordance with our above result of Eq. |3T] in one component bulk condition. Substituting 
these expressions of excess energy density and chemical potential, and those of the ideal parts into 
Eq. |36l we finally get 


— = = 1 + 4 V ??* = 1 + 477 + 477^ + .... (42) 

p 1-77 ^ 

It’s clear the first two virial coefficients from Eq.H5]is exact, which is reasonable with the approximation 
that we do not take into account the terms of order 0{rf‘){n ^ 2). For the later terms in the order of 
0{ri'^){n ^ 2) = 0{a^){m ^ 6 ) in Eq. 011 they can not be predicted well in LHS, since the terms of 
order 0{a'^){m ^ 6 ) are ignored. 

For the size modified model[^ 1^. 1^ with one component bulk fluid, Frydel and Levin have arrived 
at the following conclusion [T^ 


P 


1 ^ ^ 

-ln(l - 77 ) = 1 + V -—v" = 1 + -77 + .... 

77 7 + 1 2 


(43) 


Different from the FMT and LHS predictions, above equation only provides the first virial coefficient 
accurately. Figure [H illustrates the compressibility factor Z = ^ versus packing fraction 77 of the three 
different results discussed in Eq. |38l Eq. HU and Eq. |43l From the inset of Figure [T1 it’s notable that 
at low 77 values, the predictions between FMT and our local approach are quite close to each other, 
while the values from the size modified models SMPB/SMPNPj^ [l^, are much lower than those 
from FMT and LHS. 



Fig. 1 The equation of state of pure hard sphere fluid. 
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2.3 Numerical Method 

Finite element method (FEM) is employed in our work to solve the three different models, PNP, 
SMPNP and LHSPNP, numerically. To accelerate the calculations, the algorithms are implemented 
with the parallel adaptive finite element package PHGjd^. The molecular surface and volume meshes, 
that are necessary for FEM computation, are generated by TMSmesh[3, @1 and Tetgenjd^- 

Similar to our former work[37j, using Slotboom transformation, the modified NP equation Eg. 1321 
can be written in a symmetric form 

V • (AVq) = 0, (44) 

where Di = Die~^\ Ci = cie^' and Vi = Pqi(j) + fi Substituting the above transformed concen¬ 

tration into Poisson equation, we have 


K N 

= (45) 

i=l i=l 

It’s notable that the weak form of Eq. |U]is symmetric and linear in finite element method, and Eq. US] 
leads to a nonlinear weak form. Similar weak forms and the bilinear forms for FEM to solve the SMPNP 
have been presented explicitly in our former work[^ [^ . To deal with the nonlinearit y of the weak 
form of Eq. US] we use Newton method to linearize the system as given in previous work [01^ . Also, 
relaxation iteration is employed to obtain convergent simulation results for the coupled PDE systems. 


3 Result and Discussion 


K+ Cl- 





Fig. 2 Concentrations of K"*" and Cl at two different bulk concentrations: (1) Cb ~ O.IM (upper panel) and 
(2) Cb = 0.5M (lower panel) when the central charge of spherical cavity Q is -lOcc. 


In this section, a spherical cavity of radius lOA in ionic solution is taken as a simulation example. 
The total computational region is a sphere of radius 80A which has the same origin as the spherical 
cavity. Numerical tests are first performed in a 1:1 electrolyte solution, in which K+ and Cl“ are added, 
to study the influence of ionic size effects in various conditions by comparisons among the LHSPNP, 
SMPNP, PB models and the MC simulation. It’s notable that the PB description for continuum 
model is equivalent to the steady state PNP equations at equilibrium state[^. The diameters of 
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these two ion species required in LHSPNP, SMPNP and MC calculations are: a(K+) = 5.5lA and 
a(Cl~) = 6.37A[31|. These values are larger than pure ionic diameters because a hydration shell is 
considered around ions in solution. Figure [5] and Figure [3] illustrate the density profiles when bulk 
concentrations are O.IM (upper panels) and 0.5M (lower panels). It is worth noting that because of 
volume exclusion, the ion density in MC simulation vanishes in the region within an ion radius to the 
cavity surface. To make a clearer comparison, the curves from MC data in both figures are shifted half 
the ion diameters towards the origin. As shown in the upper panel of FigureOat low bulk concentration 
and low central charge, the simulation results among the four methods are close to each other and 
the difference between the SMPNP and LHSPNP model is small. However, when the central charge 
increases to -dOcc shown in Figure [31 the concentration of K+ from the MC simulation is considerably 
larger than those from the SMPNP and LHSPNP calculations. This seems that the SMPNP and 
LHSPNP models overestimate the volume exclusion effects in the environment of strong electric field 
and condensed ionic density. While in MC simulation, the ion correlation effects that are missed in 
SMPNP and LHSPNP can enhance the condensation of counter ions in the situation. As shown in 
both Figures [3] and [31 the LHSPNP tends to result in a slightly higher counter ion concentration near 
the sphere than the SMPNP model. 


K+ ci- 



Fig. 3 Same as in Fig|2]but for Q = —40ec. 


In order to investigate the competition among various kinds of counterions, two positive ion species, 
K+ and Ca^’*', and one negative ion species, Cl“ are added in the solution. The neutrality condition 
is applied on the boundary of the computational region. In the following, if not specified, the bulk 
concentration of K+ is set to be Cbi = O.IM. Only the bulk concentration of Ca^+, Cb 2 is varied. The bulk 
concentration of Cl“, Cbz is then determined from the neutrality condition. Three models, LHSPNP, 
SMPNP and PB models are used to simulate the spherical cavity in ionic solution, respectively. The 
concentration profiles of different ion species in the steady state are obtained from the numerical 
solutions of those models for comparison. The spherical cavity contains a central charge of -40ec and 
the diameter of Ca^+ is 4.75A[^. In the figures below, solid curve is applied for K+, and dashed curve 
for Ca2+. 

Figure [J] shows the concentration profiles of Ca^+ obtained from the three models at two differ¬ 
ent bulk values. In Figure UKa), when the bulk concentration of Ca^+ is 10“^M, the PB predictions 
overestimate the concentrations compared with those of SMPNP and LHSPNP. And the concentration 
from PB decreases quickly to smaller values than the other two predictions in a narrow region within a 
distance of 0.4A to the cavity surface. However, when Cb 2 = I0“^M, the SMPNP and LHS predictions 
are always higher than the PB, even near the spherical surface. This indicates PB model can also lead 
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Fig. 4 Concentrations of Ca^"*" at two different bulk concentrations: (a) Cb 2 ~ 10 and (b) Cb 2 = 10 '^M. 


to underestimation of concentration for counterions when its bulk concentration is small enough. This 
may be explained as follows. According to the Boltzmann distribution, the following equations hold: 

c{K+) = c{Ca^+) = 

c(C'a2+) Cb 2 

Eq. |37]can be reduced to when Cbi and Cb 2 are O.IM and respectively. The 

reduced potential u = ecP4> is about -6 under the given neutrality boundary conditions, resulting in 
the ratio ^^^ 2 +) larger than 2. Therefore, K'*' plays a leading role in neutralizing the fixed negative 
charge, while the concentration of Ca^+ may be underestimated by the PB model. Besides, it’s easy 
to see that in both subfigures of Figure 2] the predictions from PB model decrease faster than the 
others, which can be explained by Eq. 1461 This heavily exponential decrease is resulted from the quick 
enhancement of the reduced potential u in the vicinity of the spherical cavity. 


(46) 

(47) 



Fig. 5 Concentrations of at two different bulk concentrations: (a) Cb 2 = 10 and (b) Cb 2 = 10 
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Figure [5] presents the predictions of K+ at the same bulk values mentioned above for Figure E] Ac¬ 
cording to the neutrality boundary condition, the Cbs values are then 0.102M and 0.1002M, respectively. 
As is known, the counterion concentration predicted from the traditional PB model is unphysically high 
in the vicinity of the biomolecule [1, This is clearly shown in the two subhgures of Figured The 
concentration on the spherical surface can reach as high as 13.5M when Cb 2 = 10“^M. As the radical 
distance increases, the concentration decreases to the bulk value quickly. For SMPNP and LHSPNP 
models, a stratification of K’*' is observed, which is quite different from the monotonic phenomenon of 
the PB model. With the increasing of the radical distance, the concentration first increases to a highest 
value then decreases to the bulk value slowly. In SMPNP, ionic size effects are incorporated through 
adding an additional entropy term to the electrostatic energy, leading to an extra excess chemical po¬ 
tential. In LHSPNP, hard sphere repulsion described by the FMT is added to the electrostatic energy 
and then the excess chemical potential is approximated locally. Both models can capture the stratifi¬ 
cation of K’*' through adding extra terms to the chemical potential. This change originates from the 
reasonable assumption that ions in solution are regarded as finite volumes rather than point charges. 
Since Ca^^ carries one more positive elementary charge than K+, huge accumulation of Ca^’*' around 
the charged spherical cavity prevents the enhancement of K+. Nevertheless, as the radical distance 
becomes larger, the electric field strength becomes smaller and so does the concentration of Ca^’*'. 
Thus the repelling from the Ca^+ weakens and the concentration of K+ increases to a certain value 
in a finite distance. At further distance, the weakening of electric field makes the concentration of K+ 
decrease. This explains the stratification of K+ from results of SMPNP and LHSPNP models. 



Fig. 6 Concentrations of two ion species: (a) and (b) Ca^^ calculated from the original PB model. 
Four different bulk concentrations are considered: (1) PB4: Cb 2 = (2) PBS: C (,2 = 10“®M, (3) PB6: 

Cb 2 = 10“®M and (4) PB7: Cb 2 ~ 10“^M. 


Next, we vary the Cb 2 value from lO^'^M to 10“^M in which the ratio of ^ can reach as high as 10®, 
while the other conditions remain the same as aforementioned. Under these circumstances, Figures [Bl- 
[8] illustrate the concentration profiles of the two positive ion species resulted from the PB, SMPNP 
and LHSPNP models. In Figure El it’s observed that both positive ion concentrations decrease as the 
radical distance increases and that the concentration of K+ is larger than that of Ca^+ under same 
conditions. These phenomena can be explained by the Boltzmann concentrations defined in Eq. 1461 
and Eg. 1471 Similar explanations have been presented in above analysis for Figured and Figure!^ In 
addition, with the decrease of the Cb 2 , the enhancement of c(K+) is clear (Figure El^a)), while c(Ca^+) 
decreases accordingly (Figure Elb)). This also occurs in the SMPNP and LHSPNP models as shown 
in Figure [7] and Figure [5] 

In Figure[7Ka) and FigureEKa), the stratification of K+ disappears when the bulk concentration of 
Ca^’*' is smaller than 10“®M. This is a straightforward result from the decrease of Cb 2 - Though Ca^+ 
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Fig. 7 Concentrations of two ion species: (a) K"*" and (b) when ionic size effects are incorporated. Four 

different bulk concentrations are considered: (1) SMPNP4: Cb 2 = 10”“*]^, (2) SMPNP5: Cb 2 = 10“®M, (3) 
SMPNP6: C 62 = 10"®M and (4) SMPNP7: C 62 = 10"^M. 


ions are much easier to be attracted to the spherical surface than K+, the number of Ca^+ is so tiny 
that they have no significant influence on the accumulation of K’*'. Therefore, a great number of K+ 
can be attracted to the neighborhood of the spherical surface. As the radical distance increases, the 
concentration of K+ decreases gradually to the bulk values attributed to the weakening of the electric 
potential. In contrast to the stratification of K+ appeared at certain cases, the concentration of Ca^+ 
always decreases monotonically with the increasing of the radical distance, even when its bulk value is 
10“^M. This means the bulk concentration is not the essential factor for ion stratification. For SMPNP 
at given conditions, Li et al. have shown that the phenomenon depends on the ratio of ionic charge 
amount over its volumeFor the similar spherical case with mixed ion species, the concentration of 
the ion with the largest or the smallest ratio changes monotonically as the radical distance increases, 
while the other ion species’ concentrations may appear with stratification phenomenon. In our settings, 
the following relation holds: 

. q{K+) qjCl-) 

a3(Ca2+) a3{K+) a^Cl-)' ^ ’ 

As a result, the change of Ca^"*- is always monotonic, while the stratification occurs only for K'*'. In 
Figure m the numerical results from LHS models also show these properties. 

However, difference between the SMPNP and the LHSPNP models is also observed in Figure [7] and 
Figure [51 When C {,2 = lO^^M, the stratification of K+ occurs in the SMPNP model but not in the 
LHSPNP model. Furthermore, there is also a small difference between Ca^+ concentration profiles from 
SMPNP (Figure EKb)) and LHSPNP (Figure IHKb)) models. The Ca^+ concentration from LHSPNP 
simulations is smaller than that from SMPNP under the same boundary conditions. This indicates that 
the LHSPNP model captures more about the exclusion effects of Ca^+ than the SMPNP model. In 
SMPNP, when the bulk concentration Cb 2 is lO^^M, Ca^+ is still much easier to be attracted than K+. 
But for LHSPNP model, the number of Ca^"*- is not enough to prevent K+ from accumulation around 
the sphere. Therefore, in LHSPNP, the stratification of K+ does not happen under this condition. 
When Cb 2 is 10“®M or less, no stratification occurs in either SMPNP or LHSPNP modeling. 


4 Conclusion 

We have proposed a local hard sphere PNP model to account for hard sphere repulsion in three 
dimensional ionic solutions. Compared with the PNP model, the LHSPNP model contains an extra 
local excess chemical potential, derived from the expansion of the variation of the excess Helmholtz free 
energy (excess chemical potential) from FMT. This local expression avoids solving of integro-differential 
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Fig. 8 Concentrations of two ion species: (a) K"*" and (b) Ca^^ when local hard sphere excess chemical potential 
is included. Four different bulk concentrations are considered: (1) LHS4: Cb 2 = 10““^]^, (2) LHS5: Cb 2 = 10“®M, 
(3) LHS6: Cb 2 = lO'^M and (4) LHS7: Cb 2 = 


equations which requires much more computer memory and time. It is interesting to notice that one 
component of our local excess chemical potential is similar to the key part in SMPNP model. This 
indicates that in some sense, the LHSPNP model essentially contains the previous SMPB/SMPNP 
models d, [13,1131, though these two models are from different backgrounds. The closeness between 
these two models under certain conditions are also demonstrated by numerical computations in this 
work. 

Theoretical study on the equation of state of one component homogeneous fluid shows that the 
LHS model can predict exactly the first two virial coefficients, performing better than the size modified 
model which only provides the first coefficient accurately. Numerical tests for an example of a spherical 
cavity in ionic solutions indicate the LHSPNP model can avoid unphysical accumulation of counterions 
around biomolecular surface. But when the bulk concentration and the potential are high, the con¬ 
centration results from LHSPNP model are lower than those from the MC simulation. Under certain 
bulk concentrations in mixed ionic solution, we find that the concentration of one counterion species in 
LHSPNP equilibrium simulation is higher than that in PNP simulation. Furthermore, the stratification 
of counterion is observed when two different counterion species are included in the solution system. 
These phenomena from LHSPNP model are quite similar to those from SMPNP. 
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